Collapsed Gibbs Sampler with GSS Slice Sampling

1 Description

Uses a collapsed (marginalizes over a, b, amp numerically) gibbs sampler which samples from the conditionals with slice sampling. Also re-parameterizes density with (sqrt(gwidth^2+lwidth^2), gwidth-lwidth) instead of (gwidth, lwitdh) if min(gwidth/lwidth, lwidth/gwidth) > 0.05 (otherwise transformation becomes unstable). Finds global maxima in posterior using population based golden section search and uses it as the starting point for the slice sampler. About x10 times faster than the previous implementation. Using Zell Prior could be another x10.

2 Data

library(spectralmc)

#toluene data
cdata <- get_toluene_measured_data()
x <- cdata$x
y <- cdata$y

3 Algorithm

peak_spawner <- function() {
  params <- list(
          amp = rnorm(1, 1, 1), # doesnt matter
          lwidth = rnorm(1, 0.1, 0.001), # matters!!!
          gwidth = rnorm(1, 5, 0.001), # matters!!!
          pos = runif(1, 1, 1), # doesnt matter
          a = 0, # doesnt matter
          b = 0 # doesnt matter
  )
  return(params)
}

fit_kp <- 30
max_a <- 1e-2
lower <- c(-Inf, -max_a, 0)
upper <- c(Inf, max_a, Inf)
noise_sigma <- 100 # estimated by eye; intentionally to big such that tiny peaks are overnoised
signal_model <- voigt.model

res <- chains.rjmcmc_like_collapsed_slicer_gibbs(
        x,
        y,
        signal_model,
        peak_spawner,
        lower,
        upper,
        noise_sigma,
        fit_kp = fit_kp,
        add_peak_every = 10,
        iter = 20000,
        print_bar = TRUE,
        half_steps = TRUE,
        decorr_trafo = TRUE,
        save_loc = save_loc,
        checkpoint_every = 50,
)

4 Results

4.1 Fit

4.1.1 All

plt.fit.chain.interactive(x, y, res$samples, signal_model, step_width = 200)

4.1.2 Burn-in

plt.fit.chain.interactive(x, y, res$samples[seq(1, 500)], signal_model, step_width = 5)

4.1.3 Traceplots

plt.traceplot(res$samples, "pos")

plt.traceplot(res$samples, "amp")

plt.traceplot(res$samples, "lwidth")

plt.traceplot(res$samples, "gwidth")

plt.traceplot(res$samples, "a")

plt.traceplot(res$samples, "b")

4.1.4 Posterior plot

plt.posterior(res$post_vals, steps = NULL)

5 Todo

• Critical Bug amp, a, b are not sampled instead they are set to roughly their means (not really either); Easy to fix.
• Not a bug: Uses a safeguard which prevents messing up the fit if the slice sampler fails, this however destroys the gibbs sampler property and causes rejection plateaus.
• Done: nu, be transformation becomes bad (why?) if either lwidth or gwidth is really small compared to the other, should switch to nu,be only if min(lwidth/gwidth, gwidth/lwidth) > 0.05

---

1. RWTH Aachen, philipp.meissner@rwth-aachen.de↩︎